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ABSTRACT 

This work investigates the propagation of magnetized, isolated old neutron stars through the interstel- 
lar medium. We performed axisymmetric, non-relativistic magnetohydrodynamic (MHD) simulations of 
the propagation of a non-rotating star with dipole magnetic field aligned with its velocity through the 
interstellar medium (ISM). Effects of rotation will be discussed in a subsequent work. We consider two 
cases: (1) where the accretion radius Race is comparable to the magnetic standoff distance or Alfven 
radius Ra and gravitational focusing is important; and (2) where Race << Ra and the magnetized star 
interacts with the ISM as a "magnetic plow" , without significant gravitational focusing. For the first 
case simulations were done at a low Mach number Ai = i for a range of values of the magnetic field 
B*. For the second case, simulations were done for higher Mach numbers, Ai — 10, 30, and 50. In both 
cases, the magnetosphere of the star represents an obstacle for the flow, and a shock wave stands in 
front of the star. Magnetic field lines are stretched downwind from the star and form a hollow elongated 
magnetotail. Reconnection of the magnetic field is observed in the tail which may lead to acceleration 
of particles. Similar powers are estimated to be released in the bow shock wave and in the magnetotail. 
The estimated powers are, however, below present detection limits. Results of our simulations may be 
applied to other strongly magnetized stars, for example, white dwarfs and magnetic Ap stars. Future 
more sensitive observations may reveal long magnetotails of magnetized stars moving through the ISM. 

Subject headings: accretion, dipole — plasmas — magnetic fields — stars: magnetic fields — X-rays: 
stars 



1. INTRODUCTION 

There are many strongly magnetized stars moving 
through the ISM of our Galaxy. One of the most numerous 
populations is that of isolated old neutron stars (IONS) 
and old magnetars, which are not observed as radio or 
X-ray pulsars but which may still be strongly magnetized. 
There are about 1000 isolated radio pulsars observed in the 
Galaxy. The typical age of a radio pulsar is estimated as 
~ 10'' yr (e.g., Manchester & Taylor 1977). Subsequent to 
the pulsar stage, the neutron stars are still strongly mag- 
netized. Pulsar magnetic fields decay on a longer timescale 
than the lifetime of a radio pulsar. Thus the number of 
magnetized isolated old neutron stars (MIONS) should be 
larger than the number of pulsars. Total number of (mag- 
netized and non-magnetized) isolated old neutron stars is 
estimated to be 10^ — 10^. The IONS could be observed 
in the solar neighborhood owing to a low-rate accretion to 
their surface from the ISM (Ostriker, Rees, & Silk 1970); 
Schvartsman 1971; Treves & Colpi 1991, Blaes & Madau 
1993). Many of them may have strong magnetic fields, 
B ^ 10^ — 10^^ G during significant period of their evolu- 
tion ~ 10* — 10^ years (e.g. Livio et al. 1998, Treves et al. 
2000). 



Recently it has been emphasized that some neutron 
stars, termed magnetars, may have anomalously strong 
magnetic fields at their origin B ~ 10^^ — 10^® G (Dun- 
can & Thompson 1992; Thompson & Duncan 1995 (here- 
after TD95); Thompson & Duncan 1996). Magnetars pass 
through their pulsar stage much faster than classical pul- 
sars, in ~ lO'' years (TD95). Observations of soft gamma- 
ray repeaters (SGRs) and long-period pulsars in super- 
nova remnants, especially young supernova remnants (Va- 
sisht & Gotthelf 1997) support the idea that these ob- 
jects are magnetars (Kulkarni & Frail 1993; Kouveliotou 
et al 1994). The estimated birthrate of SGRs is - 10% 
of ordinary pulsars (Kulkarni & Frail 1993; Kouveliotou et 
al. 1994, 1999). Thus magnetars may constitute a non- 
negligible percentage of IONS. 

MIONS and magnetars typically move supersonically 
through the ISM and have extended magnetospheres. Two 
main regimes are possible: In the first, the Alfven radius 
Ra is much smaller than gravitational accretion radius 
Racci so that matter is gravitationally attracted by the 
star and direct accretion to a star is possible (e.g. Hoyle 
and Lyttleton 1939; Bondi 1952; Lamb, Pethick, & Pines 
1973; Bisnovatyi-Kogan & Pogorelov 1997). In the second 
regime, the magnetic standoff distance or Alfven radius Ra 
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Fig. 1. — Geometry of the MHD simulation model. The solid lines are magnetic field lines which are constant values of 
the flux function, ^'(r, z) =const. The "if values shown arc equally spaced between = 2 x 10~^ and '^max — 10^^ in 

dimensionless units discussed in §3. 



is larger than accretion radius and the magnetosphere in- 
teracts with the ISM without gravitational focusing. This 
case we term the "magnetic plow" regime (this corresponds 
to "georotator" regime in Lipunov, 1992). This is the 
regime for fast moving MIONS and magnetars. Neither of 
these regimes was investigated numerically in application 
to magnetized star propagation through the ISM. Most 
of simulations of this type were done to model the inter- 
action of the Earth's magnetosphere with the Solar wind 
(e.g., Nishida, Baker & Cowley 1998), where parameters 
of the problem were fixed by the Solar wind and Earth's 
magnetic field. 

In this paper we investigate the supersonic motion of 
magnetized stars through the ISM where a wide range of 
physical parameters is possible. We investigate the physi- 
cal process of interaction of magnetospheres with the ISM 
and estimate the possible observational consequences of 
such interaction. In §2 we estimate the important physi- 
cal parameters, and in §3 we describe the numerical model. 
In §4 we summarize results of simulations for Ra <j Race 
and for a small Mach number, A4 = 3. In §5 we discuss 
results of simulations in the "magnetic plow" regime. In 
§6 we discuss possible observational consequences of our 
results. In §7 we give a brief summary. 



2. PHYSICAL MODEL 

After the radio pulsar stage, neutron stars are still 
strongly magnetized and rotating objects. This work 
treats the motion of a non-rotating magnetized star 
through the interstellar medium. Treatment of the mo- 
tion of a rotating star through the ISM is discussed by 
Romanova et al. (2001). 

A non-magnetized star moving through the ISM cap- 
tures matter gravitationally from the accretion radius (e.g. 



Shapiro & Teukolsky 1983), 
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where V200 = w/(200 km/s) is the normalized velocity 
of the star, Cg the sound speed of the undisturbed ISM, 
and Ml. 4 = M/(1.4M0) is the normalized mass of the 
star. The mass accretion rate at high Mach numbers 
A4 = v/cs >> 1 was derived by Hoyle and Lyttleton 
(1939), 

Mhl = 47r(GAf)24 « 9.3 x lO^^M^^ ^ , (2) 

where p is the mass-density of the ISM and n = n/1 cm~'^ 
is the normalized number density. For arbitrary A4 a gen- 
eral formulae was proposed by Bondi (1952), 



Mbhl = T^RlccP"" = 47ra(GM)^ 
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where the coefficient a is of order unity (e.g., Bondi 
proposed a — 1/2, see also Rufi'ert 1994a, b; Pogorelov, 
Ohsugi, & Matsuda 2000). 

For the case of a moving magnetized star, the standoff 
distance where the inflowing ISM is stopped by the star's 
-B— fleld is referred to as the Alfven radius Ra- For a rela- 
tively weak stellar magnetic field Ra ^ Race, and in this 
limit of "gravitational accretion" denote the Alfven radius 
as Raq ■ The accretion flow becomes spherically symmetric 
inside Race, and one flnds 
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(e.g.. Lamb, et al. 1973, Lipunov 1992), where J5, is the 
magnetic field at the surface of the star of radius i?* and M 
is the accretion rate. If a magnetized star accretes matter 
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with the same rate as a non- magnetized star, M = Mbhl, 
then the Alfven radius is 
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which is about Racc/Si for the adopted reference parame- 
ters. Here, B12 = B./IO^^q ^nd = RjlO^cm. 

However, there is reason to beheve that a magnetized 
star accretes matter at a lower rate than a non-magnetized 
star for the same v, Cs, and M. Our study of spheri- 
cal Bondi accretion has shown that the magnetized star 
accretes at a lower rate than the same non-magnetized 
star (Toropin et al. 1999, hereafter T99). The magne- 
tosphere acts as an obstacle for the flow, thus decreasing 
the rate of spherical accretion compared to the Bondi rate 
Mb = 47ra(GM)2p/c3. Equations (28) and (32) of T99 
correspond to the approximate dependence 
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for R*^g/R^ in the range 1—10, where i?^^ is given by 

equation (4) with M = Mb- Thus, for a larger R'^^g, M 
is smaller and the actual Alfven radius given by equation 
(4) is larger. 

Equation (6) was deduced from simulations at small val- 
ues of Raq/R* and therefore it cannot be reliably extrap- 
olated to very large values of this ratio. Instead we can 
write in general M = ICMb where /C < 1. Then we find 
that the actual Alfven radius is Raq = R^'Ag^ '^^''- The 
two radii. Race and Raq are equal at /C w 10"'^ for our 
reference parameters. It is not known whether accretion 
can be so strongly inhibited at such small values of /C. 

Magnetars have significantly stronger magnetic fields 
than typical radio pulsars, and consequently most of them 



are in the "magnetic plow" regime. Comparison of equa- 
tions (1) and (5) shows that Race < Ra if 
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Thus, even for JC = 1 and v ~ 200 km/s, magnetars are in 
the magnetic plow regime. 

In the "magnetic plow" regime, the Alfven radius Rap 
follows from the balance of the magnetic pressure of the 
star B'^/iir = B^{R^,/R)^ against the ram pressure of the 
ISM which is pv^ for Mach numbers 3> 1. Thus 
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The magnetic field strength at this distance from the star 
is 

Ba = (47rp)i/2^ sa 9.2 x lO^^n^/^^zoo G . (9) 

At the boundary between the gravitational and magnetic 
plow regimes, equations (1), (4), and (8) coincide. We 
mention here the important influence of rotation of the 
star. Due to the rotation the magnetic field decreases with 
distance esc 1/r at large distances rather than c>c 1/r^ so 
that the Alfven radius is much larger than one described 
by equation (8) (Romanova et al. 2001). 

The velocity distribution of MIONS and magnetars is 
unknown, but it is expected to be similar to that of radio 
pulsars. Pulsars have a wide range of velocities, 10 km/s < 
V < 1500 km/s, with the peak of the distribution at w « 
175 km/s (Cordes & Chernoff 1998). Some authors give 
a smaller value, v « 100 km/s (Narayan & Ostriker 1990) 
while others give larger values, v « 250—300 km/s (Hansen 
& Phinncy 1997), v « 500 km/s (Popov et al. 2000). For 
temperatures of the ISM, T « lO'* K, the sound speed of 
gas is Cs ~ y^"fkBT/fh « 11.7 km/s, where fh « nip is 
the mean particle mass, fc^ is Boltzmann's constant, and 




Fig. 2. — Results of hydro simulations of accretion to a non-magnetized star at Mach number A4 = 3. The background 
and contours represent density. Arrows represent velocity vectors. 
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Fig. 3. — Results of simulations of accretion to a magnetized star with magnetic field B^, = 3.5 at Mach number 
= 3. Poloidal magnetic field lines and velocity vectors v are shown. The background represents density. The thick 
line represents the Alfven surface. 



7 is the usual specific heat ratio. Thus, the Mach num- 
ber of radio pulsars is in the range A4 = v/cs ~ 1 — 150 
with most pulsars having ~ 10 — 50. The accretion 
radius Race depends strongly on the velocity of the star 
which may change the ratio between Race and Ra and cor- 
respondingly the regime of accretion. For example, very 
fast MIONS with v ^ 1000 km/s have much smaller ac- 
cretion radii than slow ones, and may have Ra >> Race 
for wide range of magnetic fields. 

It is clear from the range of surface magnetic fields of 
MIONS and magnetars and the range of their velocities 
that different regimes are possible: (1) The regime of 



gravitational accretion. Race >> Ra', (2) The interme- 
diate regime. Race Ra', and (3) The "magnetic plow" 
regime. Race << Ra- In this paper, we present results 
for regimes (2) and (3), which are characterized by forma- 
tion of extended magnetotails. Regime (1) will be investi- 
gated in a future work. Below, we present numerical model 
and results of simulations, and return to discuss physical 
model further in §6, where the possible observational con- 
sequences are considered. 



3. NUMERICAL MODEL 




Fig. 4. — Same case as Figure 3, but the inner region is shown at higher resolution. Arrows show matter flux vectors 
pv. The thick line represents the Alfven surface. 
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Fig. 5. — Same case as Figure 3, but the streamlines of matter flux pv are shown. Background and dashed lines represent 
logarithm of magnetic flux which is equally spaced between logj^Q 5* — —6 and logj^Q 4* = —4. The minimum value of \1/ is 
smaller than that in Figure 3. Numbers show the value of logarithm of 



To investigate the interaction of a magnetized star with 
the ISM we use an axisymmetric resistive MHD code and 
arrange dipole so that its axis is aligned with the matter 
flow (see Figure 1). The code uses a total variation dimin- 
ishing (TVD) method (Savelyev et al. 1996; Zhukov et al. 
1993). The code was used earlier for a study of spherical 
Bondi accretion to a star with dipole magnetic field (T99). 

We used a cylindrical coordinate system (r, (f), z) with 
its origin at the star's center. The z-axis is parallel to 
the velocity of the ISM at large distances Vqo . The dipole 
magnetic moment of the star /x is parallel or antiparallel 
to the z— axis. Axisymmetry is assumed so that d/ dcf) — 
for all scalar variables. We solve for the vector potential 
A so that the magnetic field B = V x A automatically 
satisfies V • B = 0. 

The flow is described by the resistive MHD equations, 



dp 
dt 
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(13) 



The variables have their usual meanings. The equation of 
state is p = (7 — 1) pe, with 7 the specific heat ratio. In 
the simulations presented here 7 = 5/3. The equations 



incorporate Ohm's law J — cr(E + v x B/c), where a is 
the electrical conductivity. The corresponding magnetic 
diffusivity rjm = c^/iina) is taken to be a constant. 

The simulations were done inside a "cylindrical box," 
< r < Rmax)- A uniform (r, z) 
mesh was used with size iV^ x Nz- The magnetized star 
was represented by a small cylindrical box with dimen- 
sions i?* << Rmax and I << Zmax, which constitutes 
the "numerical star." In equation (11) the gravitational 
force is due to the star, = ~GM p'R./R^ . The gravita- 
tional force was smoothed for distances inside the region 
r = l^l = 0.25i?* which does not influence the computa- 
tional results outside of the numerical star. 

A point dipole magnetic fleld B = [3R (/x • R) — 
R^fi]/R^ with vector-potential A = /x x 'R./R^^ was ar- 
ranged inside numerical star at the radii: r > 0.25i?*. 
This dipole field differs from that used in T99 where a 
small but finite size "current" disk was used to produce 
the dipole field. A similar model of the field was used by 
Hayashi et al. (1996), Miller & Stone (1997), Goodson et 
al. (1997). 

The vector potential was fixed inside the numerical star 
and at its surface during the simulations. These condi- 
tions follow from the E and B boundary conditions on 
the surface of the perfectly conducting star and protect 
the magnetic field against numerical decay (T99). The 
hydrodynamic variables p, Vr, Vz, and e were fixed at the 
surface of the numerical star. These conditions are similar 
to the standard "vacuum" conditions adopted in hydro- 
dynamic simulations (e.g., Ruffert 1994a, b). However, the 
vacuum is not made too strong because of the difficulty 




Fig. 6. — Same case as Figure 5, but the inner region is shown at higher resolution. 
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Fig. 7. — Panel (a) shows the density and panel (b) the velocity variation along the z— axis. Panel (c) shows the 
energy-density variation along the z— axis, and panel (d) shows its variation with r &X z = \. Here, E^g^g—is the magnetic 
energy-density, Ekin— is the kinetic energy- density, and Egas is the thermal energy-density. For the case shown, = 3.5 
and = 3. Dotted lines on panels (a,b) correspond to hydrodynamic simulations. 



of handling low densities in MHD simulations. We discuss 
the boundary conditions on the numerical star further in 
§4.1. We tested the influence of the numerical star shape 
on our simulation results. Namely, we created an approx- 
imation of a sphere on rectangular grid and compared it 
with the cylindrical star and observed that the difference 
in the shapes has an insignificant influence on our results. 

We put the MHD equations in dimensionless form using 
the following scalings: The characteristic length is taken 
to be the Bondi radius, Rb — GM/cl^, where c^oo is 
the sound speed in the undisturbed ISM. Temperature is 
measured in units of Too, and density in units of p^c- The 
magnetic field is measured in units of the reference mag- 
netic field Bq. a reference speed is the Alfven velocity 
corresponding to a reference magnetic field Bq and den- 
sity Poo, vo = Bq/ ^Airpoo- Time is measured in units of 



to = {Zmax — Zmin)/i>oo, which is the crossing time of the 
computational region in the absence of a star. After re- 
duction to dimensionless form, the MHD equations (10) - 
(13) involve three dimensionless parameters. 





9 = 


= -7/3 


Vm = 


rira 


1 


Revo 


Rem 



where fjm is the dimensionless magnetic diffusivity, and 
Rem is the magnetic Reynolds number. Note that the 
first two parameters are dependent because of our choice 
of the length scale Rb- 

The external boundaries of the computational region 
were treated as follows. Supersonic infiow with Mach 
number A4 was specified at the upstream boundary {z = 
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Fig. 8. — Results of simulations of motion of a magnetized star with magnetic field = 14 through the ISM with Mach 
number = 3. Magnetic field lines and velocity vectors are shown. The background represents the density. The thick 
line indicates the Alfven surface. 



Zmin,0 ^ f Rmax)- At the dowustream boundary 
{z = Zmax,0 ^ r < Rmax), a "free boundary" condition 
was applied, d/dn = 0. Inflow of matter from this bound- 
ary into the computational region was forbidden. At the 
cylindrical boundary {Z„iin < z < Zmax,r = Rmax, we 
used the free boundary conditions and in some cases for- 
bid inflow to the computational region. We observed that 
the result is very similar in both cases. We checked the in- 
fluence of external boundary conditions by performing test 
simulations at different sizes of the computational region. 

The size of the computational region for most of the 
simulations was Rmax = ^Rb = 2, Zmin = -Rmax = -2, 
and Zmax = 2Rm.ax = 4. The grid Nr x Nz was 257 x 769 



in most of cases. The radius of the numerical star was 
i?* = 0.05Rb = 0.05 in most cases, but test runs were 
also done for i?* = 0.02. A number of different values of 
i?* were investigated in the purely hydrodynamic simula- 
tions (see §4.1). 

For most of our simulation runs f3 = 10^^. Therefore, 
our reference magnetic field Bq = a/SttPoo/ P is also fixed 
since poo is fixed. Thus a useful measure of the strength 
of the magnetic field is the ratio of the maximum value 
of the z— component of the field at the point r = 0.25i?* 
and z — to Bq. We denote this dimensionless field as 
i?, . We performed simulations for a range of values of i?, . 
The magnetic diffusivity was taken to be fjm = 10^® in 




Fig. 9. — Same run as Figure 8, but the inner region is shown at higher resolution. Arrows show matter flux vectors pv. 
The thick line represents the Alfven surface. 
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Fig. 10. — Same case as Figure 8, but the streamlines (solid lines) of matter flux pv are shown. The background represent 
logarithm of magnetic flux which is equally spaced between log^Q = — 6 and log^Q ^' = —4. The numbers indicate the 
logarithm of . 



most of runs, but the dependence of our solutions on fjm 
is discussed in §5.2. 

Initially, at f = the magnetic field of the star is a 
dipole field. The density and flow velocity are homoge- 
neous in the simulation region: p = p^o and v = v^o (see 
Figure 1). We investigate subsequent evolution and follow 
the evolution as long as it is needed to reach stationarity 
or quasi-stationarity. This is typically several dynamical 
time-scales. 

4. ACCRETION FOR Ra ~ Race AND M = 3 

In this section we take the Mach number to be relatively 
small, Ai — 3, so that the accretion radius Race is of the 
order of magnitude of the Alfven radius Ra- 

4.1. Hydrodynamic Simulations 

First, for reference, we did hydrodynamic simulations 
of the BHL accretion to a non-magnetized star for Mach 
number Ai = 3. We verified that the nature of the flow is 
close to that described by earlier investigators of hydrody- 
namic BHL accretion (e.g., Matsuda et al. 1991; Ruffert 
1994b). Namely, incoming matter forms a conical shock 
wave around the star. Figure 2 shows the main features 
of the flow at a late time t = 6.7to when the flow is sta- 
tionary. The opening angle of the shock wave at large 
distances from the star relative to the z— axis is predicted 
to be 6* = arcsin(l/7W), which is 6* = 19.5° for M = 3. 
Our simulations give 9 « 25°, which is larger than pre- 
dicted. However, when we performed the simulations in 
an enlarged region, Rmax = 2, 2mi„ = -2, Zmax = 4, on 
a grid NrX Nz ^ 257 x 769, we obtained 9 « 23° which 
is close to the theoretical value and similar to the value 
obtained by Ruffert (1994b). 

We calculated the accretion rate M to the numerical 
star and got a value M « 0.5Mbhl- We performed simu- 
lations using a smaller numerical star R^ = 0.02 and got 
a slightly smaller value M « OAAIbhl- This behavior 
agrees with Ruffert 's results on the dependence of M on 
numerical star size for the sizes used, i?* = 0.25Racc and 
i?* = O.lRacc (Ruffert 1994 a,b). This size dependence 
becomes negligibly small for i?* < O.lRacc- In our simu- 
lations of accretion to magnetized stars we take the larger 
value i?* = 0.05 = 0.25Racc, because it gives better reso- 
lution of the magnetic field near the star. 



We compared simulations in the region {Rmax = 2Rb, 
Zmax = ^Rb) with simulations in twice as smaller region. 
This gave only a ~ 5% decrease of accretion rate, which 
means that our standard region Rmax — 2Rb = 2 is suffi- 
ciently large to accumulate matter from the far distances, 
though simulations in smaller regions will be also suffi- 
cient. Usually, a low pressure is arranged inside the nu- 
merical star (e.g., Ruffert 1994 a,b). However, it is impos- 
sible to perform MHD simulations with very low pressure 
and density inside the numerical star. In our MHD simu- 
lations we have p = po — 1 inside the numerical star. To 
test the influence of this value, we performed simulations 
with lower densities inside numerical star, pacc = lO^^po 
and 10"'^ Po We observed, that this changed only slightly 
the accretion rate (at the level < 5%). This is connected 
with the fact that the matter density which accumulates 
around the star before accretion is much larger than po, 
so that the difference Ap = p — po is about the same for 
considered values of po- 

4.2. Accretion to a Magnetized Star 

Next, we investigated propagation of a magnetized star 
through the ISM. Simulations were performed for a num- 
ber of values of the magnetic field . We show results for 
two cases: for a relatively weak magnetic field, ~ 3.5 
(where Ra < Race)', and for a strong magnetic field, 
= 14 (where Ra > Racc)- 

Figure 3 shows the main features of the flow for a star 
with Bir = 3.5 at time t = 5tQ when the flow is stationary. 
One can see that the magnetic field of the star acts as an 
obstacle for the flow and a conical shock wave forms as 
in the hydrodynamic case with similar angle 9 as expected 
since the Mach numbers are the same. Magnetic field lines 
(with fiux values the same as in Figure 1) are slightly 
stretched by the flow, but they remain closed. Figure 4 
shows the inner region of the flow in greater detail. The 
bold line represents the Alfven surface, where the matter 
energy-density p{e + v'^/2) is equal to the magnetic energy- 
density B'^ /{Sir)- The radius of Alfven surface in z— direc- 
tion downstream at r = is Ra ~ 0.1, and in r~ direction 
at z = is Ra ~ 0.14 which is smaller than accretion ra- 
dius Race ~ 0.2. Thus, some gravitational focusing is ex- 
pected and indeed we observe density enhancement around 
the star. Figures 5 and 6 show the distribution of magnetic 
flux with the lower limit logj^Q^'mm = — 6 compared to that 
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shown in Figures 3 and 4 where logiQ^^/rnin ~ —5.3. Thus 
the apparent truncation of the magnetosphere in Figures 3 
and 4 was connected with choice of minimum plotted mag- 
netic flux. Streamlines of matter flow pv shown in Figures 
5 and 6 reveal that matter from radii r < O.lRacc, accretes 
to the star, while the rest of matter flies away. Compared 
with the non-magnetized case, the magnetic field acts as 
an obstacle for the flow, and most of the inflowing matter 
is kept away from the star. 

The matter density is strongly enhanced in the shock 
wave, but gradually decreases as it approaches the surface 
of the star where it accretes (Figure 7a). Behind the star 
(for < z < 0.4) there is also an accumulation of mat- 
ter connected with gravitational focusing by a star. Note, 
that in the case of hydrodynamic accretion, the density 
jump in front of the star (at z < 0) is much smaller, while 
behind the star (at z > 0) it is much larger. The velocity 
Vz (panel b) decreases sharply in the shock wave to small 
subsonic values, but later increases again in the polar col- 
umn. Behind the star the velocity is negative in the small 
region 0.05 < z < 0.1, where accretion occurs. The den- 
sity and velocity jumps in front of the star do not satisfy 
the standard Rankine-Hugoniot conditions, because the 
shock wave is "attached" to the magnetosphere. Matter 
cannot move freely after passage through the shock wave, 
and extra matter accumulation occurs in the shock. 

Prom panels (a) and (b) of Figure 7 it is clear that the 
rate of accretion is smaller in the case of a magnetized star 
compared to a non-magnetized star. We observed, that 
the accretion rate to a magnetized star for B^, = 3.5 is 
about 3 times smaller than that to a non-magnetized star. 
The variations of the of energy-densities along and across 
the tail (Figure 7c, d) shows that magnetic energy-density 
dominates only in a small region around the star. 



In case of a stronger magnetic field, = 14, larger mag- 
netic flux is stretched downwind (see Figure 8). Gravita- 
tional focusing is still important and density enhancement 
is observed around the star, but it is much smaller, than in 
case of weaker magnetic field = 3.5. Now, the Alfvcn 
surface has elongated structure and extends all along the 
0— axis, so that the magnetic energy-density predominates 
in the tail (see Figure 8). 

The magnetosphere around the star is larger than in 
the case _B* ~ 3.5, and the Alfvcn radius in r direc- 
tion is Ra ~ 0.26, which is larger than accretion radius 
Race = 0.2 (Figure 9). Now, all incoming matter goes 
around the magnetosphere and flies away. Streamlines of 
the flow (Figure 10) show, that no matter goes from the 
front and accretes to the back side of the star. A small 
flux of matter coming from r << Race, accretes directly 
to the upwind pole of the star. 

Figure 11a shows, that at B* = 14, compared to = 
3.5, magnetic energy-density predominates in the tail in 
the region of equatorial plane. Figure lib and also Figure 
8 and 9, show that in r— direction magnetic energy-density 
dominates in the tube with radius r ~ 0.11. 

We performed additional simulations for magnetic flelds 
strengths = 2, 7, and 11, and derived the dependence 
of the accretion rate on the magnetic field strength i?* for 
all cases. We observed that accretion rate strongly de- 
creases with increasing magnetic field (see Figure 12) as 
M ^ 5-i-3±o.o5 rpj^jg dependence reflects the fact that 
a stronger magnetic fleld of the star deflects the incoming 
ISM flow more efficiently than weaker magnetic field. 

Figure 13 shows axial variation of for different val- 
ues of B*. In all cases the magnetic field decreases very 
gradually with z: B^ ~ z~°'^^. The decrease is partially 
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Fig. 11. Panel (a) shows the energy-density variation along the tail. Panel (b) shows the variation across the tail at 
z = 1. Both panels are for B* = 14, Ai = 3. Here, Emag is the magnetic energy-density, E^in— the kinetic energy-density, 
and Egas the thermal energy- density. 



10 







1* 












DX-1 

-1.5 






■ 

V 

rx 






0.5 

logB. 





Fig. 12. — Dependence of the accretion rate to a star on the surface magnetic field B* for a star moving at Mach number 
M=3. 



connected with gradual radial expansion of the magneto- 
sphere, and partially with the reconnection of magnetic 
field lines in the tail. In the actual flow, the magnetic dif- 
fusivity is expected to be much smaller than that in the 
code. This acts to decrease the reconnection rate and in- 
crease the length of the tail. Note that the tail of Earth's 
magnetosphere extends to more than a hundred of Earth 



radii (e.g., Nishida et al. 1998). The value of the field in 
the tail is larger for larger values of B*. Even for = 3.5 
the magnetic field stretches a long distance downwind from 
the star. The Alfvcn surface in this case is small not only 
because the magnetic field is weak, but also because mat- 
ter energy-density is high. At magnetic field strengths 
< 2 — 3, however, stretching of magnetic field to the 
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Fig. 13. — Variation of the magnetic field along the tail at r = for z ;^ 0.1, =3 and different values of the magnetic 
field at the star's surface. 
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Fig. 14. — Axial distribution of density in the tail for A4 = 3 and different values of the magnetic field i?* 
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Fig. 15. — Results of simulations for = 10 and i?, = 14. 
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tail becomes suppressed. 

For > 7 {Ra ^ Race), the density in the magnetotail 
is lower than that in the incoming flow po and it decreases 
at higher B^, (see Figure 14). Thus, one can expect hollow 
magnetic tails in the case of strongly magnetized stars. 
This is connected with the fact, that magnetosphere is an 
obstacle for the flow and the tail represents the rarefaction 
region which usually forms behind an obstacle in a super- 
sonic flow (e.g. Landau & Lifshitz 1960). Furthermore, 
external matter penetrates only slowly across the magne- 
totail, because the magnetic diffusion time-scale across the 
tail is long compared with the transit time of the matter 
in the z— direction. 



5. "MAGNETIC PLOW" REGIME {Ra » Race) 



In this section we investigate interaction of magneto- 
sphere with the ISM in "magnetic plow" regime, where 
the Alfven radius Ra is much larger than accretion radius 
Race- In this limit gravitational focusing is unimportant 
and there is only direct interaction of the ISM with mag- 
netosphere of the star. For Mach numbers larger than 
about = 3 (for our set of parameters B,), the flow is 
in the "magnetic plow" regime. In this section we investi- 
gate properties of magnetotails at different Mach numbers 
(§5.1) and different magnetic diffusivities (§5.2). 

5.1. Investigation of Magnetotails at Different Mach 
Numbers M. 

In this subsection we fix the magnetic field to be = 14 
and the diffusivity 77™ = 10~^ and investigate flows at 
Mach numbers 7W = 10, = 30 and M = 50. We ob- 
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served, that at high Mach numbers M, the sharp density 
enhancement is observed in the shock cone, while the rest 
of the tail have low density (see Figures 15, 16 & 17). At 
very high Mach number A4 = 50, some kind of instabil- 
ity appears in the tail which determines its wavy behavior 
(Figure 17). This instability may be connected with high 
velocity gradient across the tail. The Alfven radius in the r 
direction and in the upwind z direction decreases at larger 
M (see also equation 8), because the flow strips deeper 
layers of magnetosphere. This also leads to higher mag- 
netic field in the tail. Reconnection is observed as in case 
of lower Mach numbers. However, the reconnection region 
is further downwind from the star at higher M. 

The axial density variations for the three cases are 
shown at Figure 18. The case with low Mach number 
= 3 is included for reference. One can see that in 
case = 10 the density in front of the star increases to 
p front = (5 ~ 6)po a-iid then decreases sharply closer to 
the surface of numerical star. At higher Mach numbers, 
the density peak is lower. The density behind the star, 
in the tail, is small ptaii ^ (10^^ — 10~^)/9o- The density 
variation across the tail at z 1 is shown at Figure 19. 
It shows that an essential part of the tail is hollow. The 
matter flux pv is much higher for higher Mach numbers 
(Figure 20), owing to higher velocities. The instability ob- 
served at = 50 may be the Kelvin-Helmholtz instability 
connected with the large gradient in the flow velocity. 

The axial magnetic field decreases slowly with distance 
behind the star, Bz ^ z 
form as in the case A4 = 
is larger at larger Mach numbers 



(Figure 21). Thus, long tails 
3. The magnetic field in the tail 



5.2. Dependence of the Flow on Magnetic Dijfusivity 

The processes of accretion and reconnection of the mag- 
netic field depend on the magnetic diffusivity fjm- The 
fact that our code explicitly includes rjm allows us to in- 
vestigate the dependence of the flows on the magnitude of 
this quantity. This is in contrast with ideal MHD codes 
where the magnetic diffusivity unavoidably arises from the 



finite numerical grid. To study the dependence on fjrm we 
fixed the magnetic field, = 14, and the Mach number, 
= 30. We made simulation runs for a range of values 
between fj„i = 10^^ and 10^^. 

We observed that at lower magnetic diffusivity, the mag- 
netic tail (the Alfven surface) is wider in the r direction. 
Figure 22 shows the variation of across the tail at z = 1. 
One can see that at small fjm = 10~^ — 10~^, regions 
with oppositely directed magnetic field are very close to 
each other, but do not reconnect. On the other hand, 
at large f/m = 10~^ — 10~^, the magnetic field is much 
smaller, because it annihilates rapidly with distance be- 
hind the star. Furthermore, note that at large fjm, matter 
is partially decoupled from magnetic field and stretching 
of magnetic field is less efficient. Figure 23 shows the de- 
pendence of axial distribution of Bz on 77^. One can see 
that at fjm > 10~^, the magnetic field decreases with z 
very rapidly. Note, that at fjm > 3 x 10"^, numerical 
diffusivity predominates, and the calculated ffows depend 
only weakly depends on fjm- 

The observed behavior is determined by the magnetic 
Reynolds number. 



RCr 



Rv Rfj 



(16) 



where tilde-quantities are our dimensionless variables. For 
example, for M. = ?> and i?, = 3.5, in the upwind re- 
gion of the flow, Rcm ~ 400 and most of the matter goes 
around the dipole and flies away or accretes to the down- 
wind pole. Matter which goes to the downwind pole has 
smaller velocity and hence smaller Rcm- Also, when mat- 
ter diffuses across the tail in the r direction owing to grav- 
itational force it has << Vz, and Rem <j 1- However, 
the timescale of the flow in the z— direction is much less 
than that in r direction, so that most of the matter flies 
away. The main conclusion of this subsection is that the 
magnetotails lengthen as the diffusivity decreases. Com- 
parison with the Earth's magnetosphere (e.g., Nishida et 
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Fig. 16. — Results of simulations for = 30 and — 14. 
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Fig. 17. — Results of simulations for = 50 and — 14. 




Fig. 18. — Axial density variations at different Mach numbers and = 14. 



al. 1998) , shows that the actual diffusivity may be smaller 
than the smallest values used in our simulations. 



6. OBSERVATIONAL CONSEQUENCES 

The question arises, is it possible to observe either a 
bow shocks or the elongated magnetotails of magnetized 
old neutron stars or magnetars? In this section we esti- 
mate the powers released and other possible observational 
features of these objects. 

6.1. Reconnection in the Tail 

Our simulations show that the magnetic field in the tail 
reconnects. This phenomenon may lead to acceleration of 
particles and possible flares in the tail. The total magnetic 



energy stored in the tail can be estimated as 

s 

Etot ^^jdz Am? [B{z)]\ (17) 



where S is the length of the tail, R{z) is the radius of the 
tail at z. The total magnetic flux in, say, the -\-z direction 
along the tail, 

$mag « B[z)AR[z)f » B R\ , (18) 

is constant in the absence of reconnection. Therefore, if 
the tail cross-section expands with distance z, then the 
magnetic field decreases as B{z) = Ba{Ra/ R{z)Y . The 
values Ra and Ba we derived earlier [see equations (8) 
and (9)]. We observed that at high Mach numbers, the 
magnetotail expands in the r direction very gradually. To 
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Fig. 19. — Radial density variations at 2; = 1 for different Mach numbers and S* = 14. 




Fig. 20. — Axial variation of the matter flux (xvz at different Mach numbers M. and = 14. 



estimate the total magnetic energy in a tail of length 5, 
we suppose that the tail does not expand, [Rtaii ~ -Ra) 
thus 

Etot ~ ^B\'kR\S ~ lO^^Bisn^^aoo 5ioo erg , (19) 

where 5*100 = ^/(lOOi?^). 

Two main physical processes determine the length of 
the magnetotail. The first is the stretching of magnetic 
field lines by the incoming matter fiow. This mechanism 
operates on the dynamical time-scale, 

S 



idyn — 



1/3 
12 



n 



-1/6,-4/3 c, 
V200 ^wos. 



(20) 



The stretched tail magnetic field has regions of oppo- 
site polarity so that the total magnetic fiux in the z di- 
rection is zero. In the axisymmetric case studied here. 



a cylindrical neutral layer forms. Magnetic field recon- 
nection/ annihilation may occur all along this layer. The 
length S of the tail is determined by the competition be- 
tween stretching and reconnection of the magnetic field. A 
nominal time-scale for reconnection across the tail is tdi / = 
R'^/rim- In view of equation (16), tdif/tdyn = Rem{R/S). 
A balance between the stretching and diffusion implies 



that this ration is of order unity. With tdif 
average power released by reconnection is 



tdyni the 



tdy 



10''B,^/^n'/'vlil ^ . (21) 



Next, we estimate the power released in an individual 
"flare," which is termed a "substorm" in the case of the 
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Earth's magnetotail. If such a flare occurs in a cyhnderical 
volume ~ ttRa^, then the energy released is 

Erec ~ ^t^Ra - lO^^Bian^/^vsoo erg . (22) 

The power of the flare, Erec = Erec /t a, depends on 
the reconnection time-scale tree = Ra/va, where va = 
Ba/^/4^ is Alfven speed. The Alfvcn speed is a function 
of density p which is uncertain. Our simulations show that 
the density in the tail is much lower than the density of in- 
coming ISM. It decreases as the magnetic field increases 
(see Figure 14). We have not been able to do simulations 
for very strong magnetic flelds such as those of magnetars. 
However, the uncertainty in utaU can be handled by look- 
ing at the extreme cases: (1) a relatively high density tail 



where ritaii — 1/cm^; and (2) a very low density tail where 
the Alfvcn velocity approaches the speed of light va ^ c. 
This density is ritaii ~ 4.4 x 10~''ni'|QQ cm~'^. 

For the case of a high matter density in the tail, we get 
Vrec ~ Va and 

tree ~ lO^Bi^/'n-VS^oo^' s , (23) 

and 

Erec ~ 10''^Bli\^/\li^o erg/s . (24) 

Note, that this power coincides with our estimate [equa- 
tion (15)] based on the dynamical time-scale. 
In the case of low density tail, we flnd 

tree ~ 7.4Bjf n-i/%2-y' s , (25) 
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Fic;. 23. Variation of the magnetic field along the tail at r = for different values of the magnetic diffusivity fjm for 
the case = 14 and M = 30. 



and the power 

Erec ~ 1.6 X lO^^'Si/^n^/^w^/^ crg/s . (26) 

Thus, the power released in individual flares is small even 
in the case of the fastest reconnection rate. The radiation 
spectrum of released energy is unknown. In view of the 
weak magnetic fields in the tail, Btaii ~ 10"** — 10~^G, 
and the possible very low densities, the energy may go 
into accelerating electrons which then radiate in the radio 
band. 

6.2. Bow Shock Radiation 

Part of the power output of a high Mach number mag- 
netized star is released in the bow show wave where the 
heated ISM behind the shock radiates. The total power 
released at the front part of the shock, r <^ Ra, is 

Eshock « IrIpv"" ~ 102in2/3^7/3^2/3 ^^^^^ ^27) 

This power is comparable to the steady power released by 
reconnection in the magnctotail. The post shock tempera- 
ture is T m.pv'^/Sk « 1.6 x 10^f2oo K which corresponds 
to X-ray band. The ISM particles excite hydrogen atoms 
which re-radiate in the optical and UV bands. Thus, one 
expects radiation from the shock wave from the optical to 
X-ray bands. 

6.3. Astrophysical Example 

In this paragraph we give the connection between 
the simulation parameters and the astrophysical quan- 
tities. The density of the ISM is taken as tIqo — 
n = 1 cm~^ and the sound speed as c^oo- Then, 
from equation (14) we obtain the reference magnetic 

field Bo ~ 0.015ni/2(c^^/30km/s)/3Z6/^ G, where 
/3_6 = /3/10~®. For example, if the dimensionless 



field is B«, then the actual magnetic field is B « 
0.015B*ni/2(c,^/30km/s)/3Zy^ G at the radius R = 
0.25i?* « 0.0125i?s = 2.6 x 10"(c^oo/30km/s)-2 cm, 
which correspond to an external region of the ac- 
tual magnctosphcre. Wc can extrapolate this field to 
smaller radii to get the magnetic field at the surface 
of the star with radius Rg = 10 km: Bg ~ 2.6 x 

10i4B*nV2(30km/s/csoo)^/3l6^^ G. 

6.4. Comparison with Earth's Magnetosphere 

There are similarities and differences between the su- 
personic solar wind interaction with the Earth's magneto- 
sphere and the interaction of the ISM with pulsars. The 
magnetization of the solar wind is important for the in- 
teraction with the Earth's magnetic field. Although not 
included in the present study, the magnetization of the 
ISM may also be important for the interactions with the 
neutron star magnetosphere. In contrast with the solar 
wind -Earth interaction, the Mach numbers of pulsars vary 
from ~ 1 to ~ 150 for the fastest pulsars (Cordes & 
Chernoff 1998). The orientation angles of magnetic axes 
9 relative to the propagation direction vary from ^ = 0° 
to ^ = 90° . If the high velocities of some pulsars are con- 
nected with initial magnetic or neutrino kicks (Lai, Cher- 
noff & Cordes 2001), then one may expect this angle to be 
closer to ^ w 0, similar to that considered in this paper. 

6.5. Observational Consequences of Long Hollow Tails 

The discussed simulations have shown that a long hol- 
low, low-density magnctotail forms behind a high Mach 
number magnetized star. This fact, and the fact that in 
the magnetic field lines are highly stretched in the tail, 
leads to possibility that particles accelerated near the star 
can preferentially propagate along the tail. This effect may 
also be important during pulsar stage. A pulsar generates 
a relativistic wind consisting of magnetic field and rela- 
tivistic particles (Goldreich & Julian 1969). The standoff 
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distance of the shock wave is determined by the total power 
generated near the Ught cyhnder (e.g., Cordes et al. 1993). 
Significant part of energy may be in magnetic field. Ex- 
panded magnetospheres of pulsars interact with the ISM 
forming elongated structures but with larger cross-sections 
compared to non-rotating stars (Romanova et al. 2001). 
Accelerated particles will propagate most easily along the 
tail of the object and may give the object an elongated 
shape. An elongated shape is observed around pulsar PSR 
2224-1-65 in the form of the Guitar Nebulae (Cordes et al. 
1993). Another elongated pulsar trail was observed in the 
X-ray band (Wang, Li & Begelman 1993). This may be 
connected with the stretching of magnetic field lines by 
the ISM. 



7. CONCLUSIONS 

Axisymmetric MHD simulations of supersonic motion of 
a star with an aligned dipole magnetic field through the 
ISM were performed for a wide range of conditions. We 
observed, that: 

1. The magnetized star acts an obstacle for the flow of 
the ISM, and a conical shock wave forms as in the hydro- 
dynamic case. 

2. Long magnetotails form behind the star, and recon- 
nection is observed in the tail. 

3. In the Ra ''^ Race regime, some matter accumulates 
around the star, but most of the matter is deflected by the 
magnetic field of the star and flies away. The accretion rate 
to the star is much smaller than that to a non-magnetized 
star. 

4. In the "magnetic plow" regime, Ra >> Race and at 
high Mach numbers M ~ 10 — 50, no matter accumulation 
is observed around the star. The density of matter in the 
tail is very low. Some matter accretes from the upwind 
pole. 



5. When Ra J> Race, the magnetic energy-density pre- 
dominates in the magnetotail. Part of this energy may 
radiate owing to reconnection processes. The power is 
however small (^ 10^^ erg/s for typical parameters for 
evolved pulsars and ~ 10^^ erg/s for magnetars), so that 
only the closest magnetars may be possibly observed. For 
tail magnetic fields B ^ 10^* — 10^® G, the tail flares or 
"substorms" may give emission in the radio band. 

6. Similar power to that released by field reconnection 
is released in the bow shock, which gives radiation in the 
band from optical to X-ray. 

7. Magnetic tails are expected to also form in case of 
propagation of pulsars through the ISM. In this case parti- 
cles accelerated by the pulsar will propagate preferentially 
along the tail to give an elongated structure. 

8. The presented simulations and estimations can also 
be applied to other magnetized stars propagating through 
the ISM, such as magnetized white dwarfs, Ap stars, and 
young stellar objects. 

9. Propagation of magnetized stars can lead to the 
appearance of ordered magnetized structures in the ISM. 
Also, these stars may give a contribution to the magnetic 
flux of the Galaxy. 

The important influence of the rotation of the star on 
the results described here has been discussed briefly by 
Romanova et al. (2001) and will be treated thoroughly in 
a forthcoming paper by our group. 
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